Classification and Characterization of the Manoor Valley’s (Lesser Himalaya) Vegetation from the Subtropical-Temperate Ecotonal Forests to the Alpine Pastures along Ecological Variables

Plant species are distributed in different types of habitats, forming different communities driven by different sets of environmental variables. Here, we assessed potential plant communities along an altitudinal gradient and their associations with different environmental drivers in the unexplored Manoor Valley (Lesser Himalaya), Pakistan. We have implemented various ecological techniques and evaluated phytosociological attributes in three randomly selected 50 m-transects within each stand (a total of 133) during different seasons for four years (2015–2018). This phytosociological exploration reported 354 plant species representing 93 different families. The results revealed that the Therophytic life form class dominated the flora, whereas Nanophyll dominated the leaf size spectra. There were a total of twelve plant communities identified, ranging from the lowest elevations to the alpine meadows and cold deserts. The maximum number of species were found in Cedrus–Pinus–Parrotiopsis community (197 species), in the middle altitudinal ranges (2292–3168 m). Our results showed that at high altitudes, species richness was reduced, whereas an increase in soil nutrients was linked to progression in vegetation indicators. We also found different clusters of species with similar habitats. Our study clearly shows how altitudinal variables can cluster different plant communities according to different microclimates. Studies such as ours are paramount to better understanding how environmental factors influence ecological and evolutionary aspects.


Introduction
The study of vegetation classification based on species co-occurrence [1,2] and its relationship to ecological variables [3] is known as phytosociology. This field has specified major

Vegetation Sampling and Herbarium work
The vegetation of the study area (Manoor Valley) was surveyed and quantified [53] during four consecutive years, from 2015 to 2018, along the environmental variables [43]. The line transect method was adopted for vegetation sampling [54][55][56][57][58]. The study area was subdivided into 133 stands (sampling plots). Each stand was replicated thrice (three transects of 50 m in each stand) [53,59]. The interval between each transect was 100 m and the interval between the stands was 200 m. The phytosociological attributes (i.e., density, frequency and their relative values, and importance value (IV)) were employed on the recorded data of each stand [4,60,61]. The species were further ranked with the highest IV and considered the representative species [19,62]. Similarly, plant communities were designated based on three dominant species [63][64][65][66]. Moreover, both attributes of the biological spectrum (life form and leaf size spectra) were recognized by the following [67]. Methods for collecting specimens, their labelling, pressing, drying, poisoning, and mounting were adopted by the following [68,69]. Their identification was achieved with the aid of Flora of Pakistan [70][71][72] and submitted to the Herbarium of Hazara University, Mansehra (Pakistan).

Ecological Variables
The slope angle, aspect and exposure were recorded at each stand using a clinometer, while the altitude, longitude, and latitude were recorded by the Global Positioning System (GPS). Two hundred grams of soil samples from three randomly selected transects within each sampling stand (0-30 cm depth) were collected [73] and mixed thoroughly to make a composite sample [74], stored in a sterile polythene bag and labeled. All the samples were submitted to the Soil and Water Testing Laboratory at the Model Farm Service Center in Mansehra, Pakistan, for analysis of various physicochemical parameters such as soil pH [75], and texture (loam, clay, silt and sand) [76], organic matter (OM%) [77], nitrogen (N) [78], potassium (K), phosphorous (P) [79], calcium carbonate (CaCO 3 ) [80][81][82], and electric conductivity (EC) [76]. Moreover, other climatic variables were measured by a small remote weather station (Kestrel 4000 weather and environmental tracker) like temperature, humidity, wind speed (WS), barometric pressure (BP), wet bulb (WB), heat index (HI), and dew point (DP) to record the data at each transect and then average values were calculated at stand level [43].

Statistical Analyses
Multivariate analysis was carried out to analyze the recorded data of species and ecological variables resulting from the field observations [49,83] to find out the relationship among them [84,85]. The recorded species and sampled stands were constrained in association to the ecological variables [86,87], which were divided into geographic, slope aspect, edaphic, and climatic variables. For the identification and classification of plant communities [53], the two-way indicator species analysis (TWINSPAN) was processed using PC-ORD version 5.0 [87][88][89]. A georeferenced map was generated with ArcGIS version 10.1 to depict the distribution pattern of plant communities.
Canonical correspondence analysis (CCA) was used to ordinate species and samples along the ecological variables [90,91] using CANOCO version 5 [92,93], and we performed a variation partitioning test (partial CCA) to evaluate how explanatory attributes (climatic, edaphic, geographic, and slope) drive the plant species distribution. First, we built the best model with the lowest number of variables (those that most explain variance), through the step function in R. Next, we also evaluated multicollinearity between variables of the final model using Variance Inflation Factor (VIF), and we removed any variable with VIF >10, one at a time. Non-multidimensional scaling ordination (NMDS) [89,94] was performed using the software R 4.0.1 [95][96][97]. NMDS was conducted to evaluate the correlation of recognized plant communities with their associated species.

Results
A total of 12 plant communities were recognized, each representing for different indicator species. Each community was associated to a set of variables, but altitude, Slope (ES), Slope (SE), Slope (SW), Slope (WN), electric conductivity (EC) and heat index were the most significant variables driving species distribution in the present study. Therophytes and Hemicryptophytes, and Nanophyll, were the most frequent type of life form and leaf size, respectively, present in the communities found. Below we discuss these results in detail.

TWINSPAN Classification
TWINSPAN, which is based partitioning reciprocal averaging ordination space, was used to classify 354 species and 133 stands. Two large different clusters, which show a high cluster heterogeneity value (Lambda = 0.814). One of these clusters had eight different communities and was formed by 93 sampling sites, while the other cluster presented four communities structured into 40 sites. Furthermore, different subdivisions observed were within these two large groups with cluster heterogeneity values of less than 0.4: a total of 12 major plant communities were recognized, from subtropical-temperate ecotonal forests (1580 m) to the alpine meadows and cold deserts (4278 m) of the Manoor Valley, Lesser Himalayas. Each community was composed of different groups of indicator species recorded at different altitudes ( Figure 1).

Vegetation Characterization of Plant Communities
In total, 12 major plant communities were established by TWINSPAN. All the twelve recognized plant communities were indicated with distinct symbols and colours. The GIS map shows the elevational layer and communities of the study area-illustrating the recognition and distribution of plant communities (

Results
A total of 12 plant communities were recognized, each representing for different indicator species. Each community was associated to a set of variables, but altitude, Slope (ES), Slope (SE), Slope (SW), Slope (WN), electric conductivity (EC) and heat index were the most significant variables driving species distribution in the present study. Therophytes and Hemicryptophytes, and Nanophyll, were the most frequent type of life form and leaf size, respectively, present in the communities found. Below we discuss these results in detail.

TWINSPAN Classification
TWINSPAN, which is based partitioning reciprocal averaging ordination space, was used to classify 354 species and 133 stands. Two large different clusters, which show a high cluster heterogeneity value (Lambda = 0.814). One of these clusters had eight different communities and was formed by 93 sampling sites, while the other cluster presented four communities structured into 40 sites. Furthermore, different subdivisions observed were within these two large groups with cluster heterogeneity values of less than 0.4: a total of 12 major plant communities were recognized, from subtropicaltemperate ecotonal forests (1580 m) to the alpine meadows and cold deserts (4278 m) of the Manoor Valley, Lesser Himalayas. Each community was composed of different groups of indicator species recorded at different altitudes ( Figure 1).

Vegetation Characterization of Plant Communities
In total, 12 major plant communities were established by TWINSPAN. All the twelve recognized plant communities were indicated with distinct symbols and colours. The GIS map shows the elevational layer and communities of the study area-illustrating the recognition and distribution of plant communities (      Figure S1). Other important variables such as altitude and windspeed (0-1.5 m/s) were found in negative association with SSI community. Nevertheless, the SSI community's species diversity was restricted by low OM (0.65-1.15%) and P (9.6 mg/kg) ( Figure 3 and Supplementary data Figure S2).

Indigofera-Parrotiopsis-Bistorta Community
The IPB community was recognized on the aspect (W-N) with a slope angle (45-80 • ) in four stands at an altitudinal range of 1789.6-1896.3 m, which has a total of 87 species associated species (Table 1) Figures S1 and S4).

Indigofera-Cedrus-Pinus Community
This plant community was recognized in 12 stands on the N-W aspect between altitudinal ranges of 1932.3-2437.8 m, which has 141 associated species (Table 1) 37.59% of all the species in the ICP community, followed by Hemicryptophytes (25.53%), Nanophanerophytes (11.35%), and geophytes (9.22%). The nanophyllous class dominated the leaf size spectra with 31.21% of species, followed by Microphyll (26.95%), and Leptophyll (17.02%). Nonetheless, aphyllous species accounted for the least number of species in the ICP community (1.42% of species, Table 1). The strongest ecological variables that significantly influenced the composition of plant species of ICP community associates were CaCO 3 (7.5 mg/kg), humidity (49.2-68.5%), and soil texture (silty loam) (Figure 3 and Supplementary data Figures S1-S4).

Abies-Picea-Juniperus Community
This community (APJ) was recorded at the middle altitudinal range (2874-3260 m) of the north-western aspect of the study area with 66 associated plant species (Table 1). The indicators of the APJ community are Abies pindrow (26.17 IV), Picea smithiana (23.77 IV), and Juniperus squamata (21.77 IV). Thymus linearis, Bistorta affinis, Bergenia stracheyi, Rheum australe, and Poa infirma are some of the herb layer's co-dominant species, while Juniperus communis and Cotoneaster microphyllus are the distinguishing species of the shrubby layer. Furthermore, Quercus incana and Pinus wallichiana are the major tree layer associates with the APJ community. The APJ community is characterised by a preference for the shade. In comparison to sub-alpine (JSJ) and alpine (SBR) communities, hill slopes get less direct sunlight. The shade effect was significantly influenced by the tree layer's larger canopy cover. Hemicryptophytes dominated the life form classes with 40.98% of species, followed by Therophytes (29.51%), Chamaephytes, and Geophytes (9.84%) each. The nanophyllous class dominated the leaf size spectra, accounting for 33.33% of plant species, followed by Microphyll (27.27%) and Leptophyll (24.24%) ( Table 1). Low K (205.4 mg/kg) and low EC (1.4 dsm −1 ) were the most significant ecological variables that played a vital role in the formation of the APJ community. Moreover, the APJ community was hosted by a clay-loamy soil texture (Supplementary data Figures S1-S4) with a low pH (Figure 3).

Sibbaldia-Bergenia-Rheum Community
This community (SBR) was identified at the higher altitudinal ranges (3199-3688 m) above the timber line at latitude (N = 34.69472-34.79333) and longitude (E = 73.60278-73.68639). The SBR community represents subalpine vegetation, having Sibbaldia procumbens (10.78 IV), Bergenia stracheyi (8.37 IV), and Rheum australe (7.75 IV) as the indicator species, for a total of 53 associated species. The herbaceous species dominated the vegetation, although some nanophanerophytes occur at comparatively lower altitudes, i.e., Juniperus squamata, J. communis, J. excelsa and Cotoneaster microphyllus. Nevertheless, other herbaceous codominants were Poa alpina, P. infirma, Bistorta affinis, and Primula hazarica. This subalpine community develops in between the timberline and alpine meadows, regardless of slope aspect, and overlaps with the alpine community (Poa-Bistorta-Primula) at most of the elevations. Hemicryptophytes dominated the life form spectra with 48% of species, followed by Therophytes (24% of species) and chamaephytes (12% of species). The Nanophyllous class dominated the leaf size spectra, accounting for 38% of all species, followed by Leptophyll (18%) and microphyll (20%) ( Table 1). The altitude and WS (2.5-5 m/s) had a significant impact on this plant community. The indicators, as well as other associated species were also found to be temperature sensitive. The soil texture hosting the SBR community was mainly clay, with the lowest K (197-216.6 mg/kg) and pH (4.8-5.8) values and maximum OM (0.98-2.28%) concentration ( Figure 3). Moreover, the SBR community was found to be negatively associated with ecological variables such as humidity, CaCO 3 , HI, WB, BP, and slope angle (Supplementary data Figures S1, S3 and S4).

Poa-Bistorta-Primula Community
The Poa-Bistorta-Primula Community (PBP) was recognized as the highest altitudinal (3724-4278 m) alpine and cold desert plant community recorded at a latitude (N = 34.69306-34.83861) and a longitude (E = 73.60750-73.69444). Poa alpina (17.32 IVI), Bistorta affinis (15.03 IV), and Primula rosea (9.07 IV) are the indicator species of PBP community. Other co-dominant species are Rheum australe, Bergenia stracheyi and Androsace hazarica. This plant community included a total of 39 species. However, tree and shrub layers (Phanerophytes and Nanophanerophytes) were entirely absent from these alpine meadows. Moreover, this alpine community (PBP) has a low species richness as compared to other plant communities ( Figure 2 and Table 1). Extremely low temperatures are a hallmark of the growth period due to high elevation. Xeric conditions compounded such harsh environments, and a relatively short growth season was recorded from July to September. Hemicryptophytes contributed 53.85% of the species, followed by Therophytes (23.08%), Chamaephytes (12.82%), and Geophytes (10.26%). Microphyll and Nanophyll classes dominated the leaf size spectra with 28.95% of species each, followed by Leptophyll (21.05%) ( Table 1). Higher altitude and WS (3.5-8 m/s), as well as low temperature had a strong impact on the indicators and other associates of PBP community. The soil texture hosting this community was sandy in nature ( Figure 3), with the lowest K (196-206.9 mg/kg) and pH (4.9-5.6) and maximum OM (1.15-2.64%) concentration. Furthermore, the PBP community was found to be negatively associated with ecological variables such as humidity, CaCO 3 , HI, WB, BP, and slope angle (Supplementary data Figures S1-S4).

Non-Metric Multidimensional Scaling (NMDS)
All the data based on 354 plant species in 133 sampled stands were categorized into 12 plant communities using NMDS. Plant communities that are close together or on the same axis have a positive correlation, whereas communities that are far apart or on different axes have a negative correlation. The Poa-Bistorta-Primula, Sibbaldia-Bergenia-Rheum, and Juniperus-Sibbaldia-Juniperus communities, for example, had a positive association with one another but negatively correlated with the Indigofera-Juglans-Isodon, and Indigofera-Cedrus-Pinus communities (Figure 4). All these relationships could be attributed to patterns in the variables of their host environment. For example, the former plant communities were found at elevations of 3724-4278 m, 3199-3688 m, and 3250-3644 m, respectively, while the latter plant communities were found at lower elevations (1597-2456 m and 1932.3-2437.8 m), respectively. Moreover, four stands (S38, S39, S40, S41 and S42) are correlated with each other and shaped the SSI community. The PBP community identified in 13 stands can be seen far apart from the other plant communities. Only the three most representative plant species for each community are plotted. attributed to patterns in the variables of their host environment. For example, the former plant communities were found at elevations of 3724-4278 m, 3199-3688 m, and 3250-3644 m, respectively, while the latter plant communities were found at lower elevations (1597-2456 m and 1932.3-2437.8 m), respectively. Moreover, four stands (S38, S39, S40, S41 and S42) are correlated with each other and shaped the SSI community. The PBP community identified in 13 stands can be seen far apart from the other plant communities. Only the three most representative plant species for each community are plotted. Figure 4. NMDS ordination is based on plant species that are in association with sampled stands and grouped into communities Plant communities that are close together or on the same axis and have a positive correlation. Only the three most representative plant species for each community are plotted. Some species are present in different communities and due to that, the total number of species are not 36.

Canonical Correspondence Analysis (CCA)
The CCA and variation partitioning tests showed that the total inertia results of CCA was 8.626, where our final variables (Altitude, Slope ES, Slope SE, Slope SW, Slope WN, EC, heat index) together explained 22.6% of variation. CCA model was significant (χ 2 = 1.951; pseudo-F value = 4.529; p < 0.001). For the 8 explanatory variables, we tested simple term effects. Simple term effects showed that all variables were significant (χ 2 range = 0.104-0.785; pseudo-F value[range] = 1.93-14.59; p < 0.006; Table 1). Finally, the two first axis were also highly significant (p < 0.001). Table 2

Canonical Correspondence Analysis (CCA)
The CCA and variation partitioning tests showed that the total inertia results of CCA was 8.626, where our final variables (Altitude, Slope ES, Slope SE, Slope SW, Slope WN, EC, heat index) together explained 22.6% of variation. CCA model was significant (χ 2 = 1.951; pseudo-F value = 4.529; p < 0.001). For the 8 explanatory variables, we tested simple term effects. Simple term effects showed that all variables were significant (χ 2 range = 0.104-0.785; pseudo-F value [range] = 1.93-14.59; p < 0.006; Table 1). Finally, the two first axis were also highly significant (p < 0.001). Table 2

Discussion
Plant species are distributed in diverse types of habitats, forming different communities driven by different sets of environmental variables [6]. The aspect stimulates habitat diversification and promotes micro-environmental variation in the vegetation structure [36,37]. As a result, the composition of different units is observed as a reflection of changing habitat settings along environmental variables [38]. Here, we used different multivariate approaches to assess potential plant communities along an altitudinal gradient and their association with different environmental drivers. Our study documented 354 plant species belonging to 93 families. The current research area is located at the elevations ranging from 1580 m to 4278 m, with varying environmental conditions that are reflected in a rich and diverse flora. Our results showed that at high altitudes, species richness was reduced, whereas an increase in soil nutrients was linked to progression in vegetation indicators. We also found different clusters of species with similar habitats. Our study clearly shows how altitudinal variables can cluster different plant communities according to different microclimates.
In the current vegetational sampling of a remote valley (Manoor Valley, Himalaya), 12 major plant communities were established by TWINSPAN from the lower ranges to the alpine meadows. The CPP communities (197 species) in the middle altitudinal habitats (2292-3168 m) have the most plant species. Ordination methods have commonly been used to show species distribution and community structure along ecological variables [98,99]. Similarly, a researcher investigated the vegetation of the western Himalayas and identified five distinct communities, the most abundant of which were found on north-facing slopes at middle altitudes, where the moisture levels were highest [6]. Thirteen major groups were identified in the vegetation of Kammanassie areas using the TWINSPAN classification [100]. These results, along with ours, show evidence that elevational variables are suitable places to evaluate how changes in environmental variables drive plant community structure and diversity. In addition, these results also show that multivariate approaches are powerful tools for community analysis [43] and can be considered in new ecological studies as statistical methods.
The vegetation in the upper altitudinal ranges includes Pinus wallichiana, Abies pindrow, Indigofera heterantha and Viburnum grandiflorum, which are the representatives of moist temperate forests. These plant species are the temperate zone representatives [30,39,41,105,106]. These plant associations were shaped by the impact of various environmental gradients. Ecosystems respond to numerous simultaneous changes in the environment as these variations differ the diversity and distribution of communities [107,108]. Vegetation in distributions more closely resembles the changes in soil characteristics [109][110][111]. Our results revealed that soil characteristics such as EC, pH, soil texture, OM, K and P had a great impact on plant community distribution and association. Soil variables, altitude, latitude, slope aspect and angle also had a strong influence on species richness, as previously reported by [112]. Dissimilar plant communities were described as those with less than 65% similarity [113,114]. The communities' similarities were due to shrubs, trees, and perennial plants, while the communities' dissimilarities were due to Therophytes. The maximum similarity index was noted between SBR and JSJ communities (55.53%), followed by JSJ and APJ communities with 39.71% of similarity. The highest similarity between communities may be due to similar environmental conditions [30], which leads to changes in the species' habitat. The highest dissimilarity was observed between PBP and IPB communities, JSJ and SSI communities (99.94% each), followed by IJI and PBP communities (99.92%), SSI and PBP communities (99.91%). These results follow the findings of [6,53]. Maximum dissimilarity between communities might be due to wide altitudinal variation among communities [30,115], which represents the presence of different set of species adapted to different set of climatic variables [112,[116][117][118].

Conclusions
To the best of our knowledge, this is the only valley within the Himalayas of Pakistan that has never been explored before, due to its harsh terrain and geographical location. The current study revealed that the sampled area has rich species diversity. The study provides the first ever detailed insights into the spatial distribution and vegetation mapping in response to environmental variables in the study area. The flora of the Manoor Valley consists of 354 plant species belonging to 93 families, distributed into a total of 12 major plant communities, from the lowest altitude to the alpine zones. The Cedrus-Pinus-Parrotiopsis community resided at the middle altitudinal ranges (2292-3168 m) was recorded with highest number of associates (197 species). Our study clearly shows how altitudinal variables can cluster different plant communities according to different microclimates, which can be a proxy for future studies evaluating the impacts of climate change on plant communities. Studies such as ours are paramount to better understand how environmental factors influence ecological and evolutionary aspects.  Figure S2: Canonical correspondence analysis: (a). the contour plot shows the count of species at each axis, (b). Distribution of plant species along the slope aspects. Figure S3: Canonical correspondence analysis: (a) Distribution of stands along the edaphic variables, (b) and (c) Van Dobben circles show the correlation of species in association to edaphic variables and among them, i.e., red circle shows the positive and blue indicates the negative correlation. Figure S4: Canonical correspondence analysis: (a) Association of sampling sites along the climatic variables, (b). Association of sampling sites along the slope aspects.